#code for DTU analysis can be found in docs DTU.Rmd. Only follow up analyses are showed in html.
rm(sumExp)
# summarise
# empirical FDR < 0.05
map_df(all_res, ~{group_by(.x, FDR = empirical_FDR < 0.05) %>% tally() %>% filter(FDR == TRUE) }, .id = "contrast")# regular FDR < 0.05
map_df(all_res, ~{group_by(.x, FDR = regular_FDR < 0.05) %>% tally() %>% filter(FDR == TRUE) }, .id = "contrast")dtu_genes <- map(all_res, ~{ .x %>% filter(empirical_FDR < 0.05) %>% pull(gene) })
#UpSetR::upset(fromList(dtu_genes),nintersects = NA)
UpSetR::upset(fromList(dtu_genes),nintersects = NA,nsets = 7)sig_res <- map_df(all_res, ~{.x %>% filter(empirical_FDR < 0.05, abs(estimates) > 1)}, .id = "contrast")
age_res <- filter(all_res$age, empirical_FDR < 0.05)
stimulation_res <- map_df(all_res, ~{.x %>% filter(empirical_FDR < 0.05, abs(estimates) > 1)}, .id = "contrast" ) %>% filter(contrast != "age")all_res %>%
map_df( ~{.x %>% filter(pval < 0.05)}, .id = "contrast") %>%
filter(!contrast %in% "age") %>%
ggplot(aes( x = estimates, y = -log10(pval), colour = empirical_FDR < 0.05 )) +
geom_point(size = 1) + scale_colour_manual(values = c("black", "red")) +
facet_wrap(~contrast, scales = "free") +
xlim(-4,4) +
theme_bw()## Warning: Removed 4 rows containing missing values (geom_point).
# age only
all_res %>%
map_df( ~{.x %>% filter(pval < 0.05)}, .id = "contrast") %>%
filter(contrast == "age") %>%
ggplot(aes( x = estimates, y = -log10(pval), colour = empirical_FDR < 0.5 )) +
geom_point(size = 1) + scale_colour_manual(values = c("black", "red")) +
facet_wrap(~contrast, scales = "free") +
xlim(-0.1,0.1) +
theme_bw()## Warning: Removed 6727 rows containing missing values (geom_point).
markers = "~/Documents/MiGASti/docs/2nd_pass/TYROBP.xlsx"
markers = read_excel(markers, col_names = TRUE)
markers = as.data.frame(markers)
isoform_id <- rownames(tx_ratio)
copy <- cbind (isoform_id, tx_ratio)
ratio <- merge(copy, markers, by = "isoform_id")
sample <- colnames(ratio)
rownames(ratio) = ratio$isoform_id
metadata <- read.table("~/Documents/MiGASti/Databases/metadata.txt")
#set rownames to Sample
row.names(metadata) <- metadata$Sample
setwd("~/Documents/MiGASti/Databases")
exclude <- read.table("samples2remove.txt")
exclude <- exclude$x
ratio = ratio[, !colnames(ratio) %in% exclude]
#Excludes the samples from filters.
#dim(genes_tpm_filt)
metadata_filt = metadata[ !(rownames(metadata) %in% exclude), ]
length(metadata_filt)## [1] 38
samples_LPS = metadata_filt$Sample[metadata_filt$Stimulation %in% "LPS"]
samples_unstim = metadata_filt$Sample[metadata_filt$Stimulation %in% "unstim"]
metadata4deg = metadata_filt[rownames(metadata_filt) %in% c(samples_LPS,samples_unstim) ,]
ratio_LPS_unstim = as.data.frame(ratio[, colnames(ratio) %in% c(samples_LPS, samples_unstim)])
ratio_m <- melt(as.matrix(ratio_LPS_unstim))## Warning in melt(as.matrix(ratio_LPS_unstim)): The melt generic in data.table
## has been passed a matrix and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(as.matrix(ratio_LPS_unstim)). In the next version, this warning
## will become an error.
gene2check_charac = merge(ratio_m, metadata4deg, by.x = "Var2", by.y = "Sample")
# show direction of effect for male versus female
ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) +
geom_boxplot(notch = F, na.rm = T) +
geom_jitter() +
scale_y_continuous(limits = c(0,1)) +
theme_bw() + facet_wrap(~Var1, scales = "free_y") +
ggpubr::stat_compare_means(label = "p.format", label.x.npc = "centre", method = "wilcox.test")## Warning: Removed 89 rows containing missing values (geom_point).
p <- ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) +
geom_violin(width = 0.5) +
geom_jitter(width = 0.25, aes(colour = Stimulation), size = 1) +
geom_boxplot(fill = NA, outlier.colour = NA, width = 0.25) +
facet_wrap(~Var1) + theme_classic()
p IFNy <- all_res$IFNy
markers = "~/Documents/MiGASti/docs/2nd_pass/WARS.xlsx"
markers = read_excel(markers, col_names = TRUE)
markers = as.data.frame(markers)
isoform_id <- rownames(tx_ratio)
copy <- cbind (isoform_id, tx_ratio)
ratio <- merge(copy, markers, by = "isoform_id")
sample <- colnames(ratio)
rownames(ratio) = ratio$isoform_id
metadata <- read.table("~/Documents/MiGASti/Databases/metadata.txt")
#set rownames to Sample
row.names(metadata) <- metadata$Sample
setwd("~/Documents/MiGASti/Databases")
exclude <- read.table("samples2remove.txt")
exclude <- exclude$x
ratio = ratio[, !colnames(ratio) %in% exclude]
#Excludes the samples from filters.
#dim(genes_tpm_filt)
metadata_filt = metadata[ !(rownames(metadata) %in% exclude), ]
length(metadata_filt)## [1] 38
samples_IFNy = metadata_filt$Sample[metadata_filt$Stimulation %in% "IFNy"]
samples_unstim = metadata_filt$Sample[metadata_filt$Stimulation %in% "unstim"]
metadata4deg = metadata_filt[rownames(metadata_filt) %in% c(samples_IFNy,samples_unstim) ,]
ratio_IFNy_unstim = as.data.frame(ratio[, colnames(ratio) %in% c(samples_IFNy, samples_unstim)])
ratio_m <- melt(as.matrix(ratio_IFNy_unstim))## Warning in melt(as.matrix(ratio_IFNy_unstim)): The melt generic in data.table
## has been passed a matrix and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(as.matrix(ratio_IFNy_unstim)). In the next version, this warning
## will become an error.
gene2check_charac = merge(ratio_m, metadata4deg, by.x = "Var2", by.y = "Sample")
# show direction of effect for male versus female
ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) +
geom_boxplot(notch = F, na.rm = T) +
geom_jitter() +
scale_y_continuous(limits = c(0,1)) +
theme_bw() + facet_wrap(~Var1, scales = "free_y") +
ggpubr::stat_compare_means(label = "p.format", label.x.npc = "centre", method = "wilcox.test")## Warning: Removed 9 rows containing non-finite values (stat_compare_means).
## Warning: Removed 243 rows containing missing values (geom_point).
p <- ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) +
geom_violin(width = 0.5) +
geom_jitter(width = 0.25, aes(colour = Stimulation), size = 1) +
geom_boxplot(fill = NA, outlier.colour = NA, width = 0.25) +
facet_wrap(~Var1) + theme_classic()
p ## Warning: Removed 9 rows containing non-finite values (stat_ydensity).
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).
## Warning: Removed 9 rows containing missing values (geom_point).
IFNy <- all_res$IFNy
markers = "~/Documents/MiGASti/docs/2nd_pass/CIITA.xlsx"
markers = read_excel(markers, col_names = TRUE)
markers = as.data.frame(markers)
isoform_id <- rownames(tx_ratio)
copy <- cbind (isoform_id, tx_ratio)
ratio <- merge(copy, markers, by = "isoform_id")
sample <- colnames(ratio)
rownames(ratio) = ratio$isoform_id
metadata <- read.table("~/Documents/MiGASti/Databases/metadata.txt")
#set rownames to Sample
row.names(metadata) <- metadata$Sample
setwd("~/Documents/MiGASti/Databases")
exclude <- read.table("samples2remove.txt")
exclude <- exclude$x
ratio = ratio[, !colnames(ratio) %in% exclude]
#Excludes the samples from filters.
#dim(genes_tpm_filt)
metadata_filt = metadata[ !(rownames(metadata) %in% exclude), ]
length(metadata_filt)## [1] 38
samples_IFNy = metadata_filt$Sample[metadata_filt$Stimulation %in% "IFNy"]
samples_unstim = metadata_filt$Sample[metadata_filt$Stimulation %in% "unstim"]
metadata4deg = metadata_filt[rownames(metadata_filt) %in% c(samples_IFNy,samples_unstim) ,]
ratio_IFNy_unstim = as.data.frame(ratio[, colnames(ratio) %in% c(samples_IFNy, samples_unstim)])
ratio_m <- melt(as.matrix(ratio_IFNy_unstim))## Warning in melt(as.matrix(ratio_IFNy_unstim)): The melt generic in data.table
## has been passed a matrix and will attempt to redirect to the relevant reshape2
## method; please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(as.matrix(ratio_IFNy_unstim)). In the next version, this warning
## will become an error.
gene2check_charac = merge(ratio_m, metadata4deg, by.x = "Var2", by.y = "Sample")
# show direction of effect for male versus female
ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) +
geom_boxplot(notch = F, na.rm = T) +
geom_jitter() +
scale_y_continuous(limits = c(0,1)) +
theme_bw() + facet_wrap(~Var1, scales = "free_y") +
ggpubr::stat_compare_means(label = "p.format", label.x.npc = "centre", method = "wilcox.test")## Warning: Removed 4 rows containing non-finite values (stat_compare_means).
## Warning: Removed 19 rows containing missing values (geom_point).
p <- ggplot(gene2check_charac, aes(x = Stimulation, y = value, fill = Stimulation)) +
geom_violin(width = 0.5) +
geom_jitter(width = 0.25, aes(colour = Stimulation), size = 1) +
geom_boxplot(fill = NA, outlier.colour = NA, width = 0.25) +
facet_wrap(~Var1) + theme_classic()
p ## Warning: Removed 4 rows containing non-finite values (stat_ydensity).
## Warning: Removed 4 rows containing non-finite values (stat_boxplot).
## Warning: Removed 4 rows containing missing values (geom_point).
length(which(all_res$LPS$empirical_FDR < 0.05))## [1] 79
length(which(all_res$LPS$empirical_FDR < 0.10))## [1] 119
length(which(all_res$LPS$empirical_FDR < 0.15))## [1] 168
LPS_t <- subset(all_res$LPS, empirical_FDR < 0.05)
res_name_LPS_top = LPS_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_LPS_top)length(which(all_res$IFNy$empirical_FDR < 0.05))## [1] 65
length(which(all_res$IFNy$empirical_FDR < 0.10))## [1] 96
length(which(all_res$IFNy$empirical_FDR < 0.15))## [1] 111
IFNy_t <- subset(all_res$IFNy, empirical_FDR < 0.05)
res_name_IFNy_top = IFNy_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_IFNy_top)length(which(all_res$R848$empirical_FDR < 0.05))## [1] 59
length(which(all_res$R848$empirical_FDR < 0.10))## [1] 80
length(which(all_res$R848$empirical_FDR < 0.15))## [1] 115
R848_t <- subset(all_res$R848, empirical_FDR < 0.05)
res_name_R848_top = R848_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_R848_top)length(which(all_res$TNFa$empirical_FDR < 0.05))## [1] 17
length(which(all_res$TNFa$empirical_FDR < 0.10))## [1] 29
length(which(all_res$TNFa$empirical_FDR < 0.15))## [1] 83
TNFa_t <- subset(all_res$TNFa, empirical_FDR < 0.05)
res_name_TNFa_top = TNFa_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_TNFa_top)length(which(all_res$DEX$empirical_FDR < 0.05))## [1] 2
length(which(all_res$DEX$empirical_FDR < 0.10))## [1] 3
length(which(all_res$DEX$empirical_FDR < 0.15))## [1] 5
DEX_t <- subset(all_res$DEX, empirical_FDR < 0.05)
res_name_DEX_top = DEX_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_DEX_top)length(which(all_res$IL4$empirical_FDR < 0.05))## [1] 48
length(which(all_res$IL4$empirical_FDR < 0.10))## [1] 81
length(which(all_res$IL4$empirical_FDR < 0.15))## [1] 147
IL4_t <- subset(all_res$IL4, empirical_FDR < 0.05)
res_name_IL4_top = IL4_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_IL4_top)length(which(all_res$ATP$empirical_FDR < 0.05))## [1] 509
length(which(all_res$ATP$empirical_FDR < 0.10))## [1] 734
length(which(all_res$ATP$empirical_FDR < 0.15))## [1] 975
ATP_t <- subset(all_res$ATP, empirical_FDR < 0.05)
res_name_ATP_top = ATP_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_ATP_top)length(which(all_res$age$empirical_FDR < 0.05))## [1] 34
length(which(all_res$age$empirical_FDR < 0.10))## [1] 108
length(which(all_res$age$empirical_FDR < 0.15))## [1] 199
age_t <- subset(all_res$age, empirical_FDR < 0.05)
res_name_age_top = age_t[, c("isoform_id", "gene_id", "gene", "estimates", "se", "df", "pval", "regular_FDR", "empirical_pval", "empirical_FDR")]
createDT(res_name_age_top)